Scaling and localization lengths of a topologically disordered system 



Jacob J. Krich 1 ' 2 and Alan Aspuru-Guzik 2 

1 Harvard University Center for the Environment, Cambridge, MA 02138 
2 Department of Chemistry and Chemical Biology, Harvard University, Cambridge, MA 02138 

(Dated: January 19, 2011) 

We consider a noninteracting disordered system designed to model particle diffusion, relaxation 
in glasses, and impurity bands of semiconductors. Disorder originates in the random spatial dis- 
tribution of sites. We find strong numerical evidence that this model displays the same universal 
behavior as the standard Anderson model. We use finite-size-scaling to find the localization length as 
a function of energy and density, including localized states away from the derealization transition. 
Results at many energies all fit onto the same universal scaling curve. 
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The disorder-induced transition from extended to lo- 
calized states in non-interacting quantum systems has 
been a rich source physics insight for over fifty years pp. 
It is relevant for a broad range of transport properties 
[J], glass formation [3J, conductivity of composites [I], 
random walks [S], as well as for non-radiative recombi- 
nation in intermediate-band photovoltaics . Scanning- 
tunneling and Bose-Einstein condensate experiments are 
increasingly able to probe localization properties directly 
[7] . The delocalization transition is usually studied using 
the standard Anderson model, which considers a non- 
interacting tight-binding lattice with uniform nearest- 
neighbor coupling and random on-site energies pQ. For 
systems in which the disorder originates in the random 
configuration of the sites rather than, e.g., random lo- 
cal fields, it is better to consider the so-called topologi- 
cally disordered or Lifshitz model, in which sites are dis- 
tributed randomly in space with no on-site energies and 
all pairs of sites are connected by hopping terms with 
amplitude exponentially decaying with the distance be- 
tween them [5J. As the density of sites increases, a lo- 
calization/delocalization transition occurs, just as in the 
standard Anderson model. The density of states [3J |5]- 
[T3] and localization properties (3j [TOj [14] of this model 
have received much attention. Here we obtain the lo- 
calization lengths of this model quantitatively as a func- 
tion of wavefunction energy and site density. We adapt 
the finite-size scaling method, which has been success- 
fully applied to the standard Anderson model [ToTll9j . 
and give strong numerical evidence that the topologically 
disordered model displays the same universal behavior. 

Model We consider particles confined to identical lat- 
tice sites distributed randomly in space, with density p' . 
Particles can hop between sites with an exponentially de- 
caying hopping coefficient. The Hamiltonian is 

H = -Y E^ K l n )( ffl l (!) 

where \n) are the site wavef unctions, r nm is the distance 
between sites n and m, a* B is an effective Bohr radius 



giving the decay of the wavef unctions, Vq is an overall en- 
ergy scale, and the sum runs over all pairs of sites. This 
model is called topologically disordered because there is 
no ordered structure describing which lattice sites are 
strongly coupled to each other. This model can describe 
the Hamiltonian of impurities with hydrogenic wavefunc- 
tions in a semiconductor, where the hopping originates 
in overlaps of their effective atomic wavefunctions. The 
dimensionless density is p = p'a*^. Experimentally, such 
systems are found to have a metal-insulator transition 
at p 1 / 3 ~ 0.26 20 . If the diagonal components of H 
are set such that each column sums to zero, the model 
can describe glass dynamics EI] and continuous time 
quantum walks [5], but many properties are similar. 

The density of states (DOS) is shown in the inset of 
Fig. [lb. The DOS is well-understood in the low-density 
[3J and high-density, low energy [HI H2] limits. 

Localization Lengths from Quasi-ID Scaling The local- 
ization length A of a localized eigenstate is determined by 
the asymptotic decay of the wavefunction tp ~ e ~' T '~' ro ' / ' A , 
where ro is some location of high wavefunction ampli- 
tude. Our goal is to find the localization length as a 
function of p and E. Finite-size scaling techniques allow 
us to study computationally tractable small systems and 
systematically extrapolate to results for the true infinite 
system. We assume, and find strong numerical evidence 
to confirm, that the topologically disordered system is 
controlled by the same fixed point as the standard An- 
derson model. The critical exponent and other critical 
parameters of this fixed point have been well-determined 
by previous studies [TTJ [TSJ [3T] , and we will use these re- 
sults to aid our study of the localization lengths. While 
studies interested in the critical exponents have focused 
on systems close to the localization transition, we are 
interested in localization lengths across the range of en- 
ergies and densities relevant to experiments, which in- 
cludes not only values close to the transition but also 
ones in the strongly localized regime. As we move away 
from the critical point, corrections to scaling should be- 
come more important, and we test our results' sensitivity 
to these corrections. 
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FIG. 1. Scaling function for data with w > 25, including one 
correction to scaling 77 with y = 2.7. The thin solid line is 
the fit. Statistical errors in A m are smaller than point sizes. 
There are 5829 data points and 893 fit parameters. The nor- 
malized x 2 statistic per degree of freedom d is 7.5. 70% of the 
data points are within error of the fit. In (a), the dashed line 
shows the asymptote a(x) — vsinh~ 1 (x/2), showing good 
agreement in the strongly localized regime. In (b), the rel- 
evant function ho(x) minus the asymptotic a(x) is plotted. 
Dashed lines show the three Gaussians from the fit of S. (in- 
set to a) A wire of width w is constructed by sequentially 
adding slices (shown at right) to an existing wire. Each slice 
has length Lq > L c and has randomly distributed particles 
of the chosen density, (inset to b) Density of states for three 
site densities, calculated with 1500 sites and 100 realizations 
of disorder. Curves offset for clarity. At high density, the 
DOS becomes asymmetric. 



The quasi-lD scaling method was introduced for the 
standard Anderson model by MacKinnon and Kramer, 
who considered several wires of varying widths w and 
long lengths [15] . We adapt the recursive Green function 
technique to find the localization length A„ for equivalent 
wires in our system, see inset to Fig.[l^, [T^IIS]- The key 
idea is to divide a long wire into slices, where the sites 
in each slice are directly coupled only to each other and 
to sites in the immediately adjacent slices. This is not 
strictly possible for the system of Eq. [I] as all sites are 
directly coupled. We take, however, a cutoff L c and set 
exp(— r/a* B ) — > for r > L c . We start with a single 
slice of width w, length L > L c and periodic boundary 
conditions in two dimensions. We recursively add slices 
adjacent to the wire and find the portion of the Green 
function that connects any site in slice 1 to any site in 
slice N at fixed energy E, We then use the standard 
method to determine the localization length X w of the 



long wire at energy E [TH [THl HH] . The width w of each 
wire must be large enough that there is a sufficient prob- 
ability that each slice will have sites within L c of sites 
in the neighboring slices, otherwise the wire is discon- 
nected [35]. Due to the varying number of particles in 
each slice and the varying off-diagonal matrix elements, 
the transfer matrix method (reviewed in Ref. I19j) cannot 
be used. 

The statistical error in A^, can be estimated by assum- 
ing that the estimate of X w after N slices is the mean of 
N independent and identically distributed samples cho- 
sen from a normal distribution, so the sampling error 
goes down as N -1 / 2 [19]. This is rigorously proved for 
the standard Anderson model with Oseledec's theorem 
[19"] , and we assume that similar statistics hold here. We 
choose a maximum error of 1%, which requires 10 4 to 10 6 
slices, depending on the parameters. 

After determining A^ for a range of widths w, MacKin- 
non and Kramer (and many subsequent authors) use the 
one-parameter scaling theory |23| to extrapolate to infi- 
nite size systems. The scaling can be expressed in terms 
of the relevant variable <j>, with <fi greater (less) than zero 
for extended (localized) states. Then the dimensionless 
quantity A K1 = X w /w is a universal function of 4> and w, 
with critical exponent v, |17j 

A w = m(p,E)w 1 ^], (2) 

The correlation length is £ = |</>|~". 

The hypothesized universality of the metal-insulator 
transition with time reversal symmetry has been con- 
firmed in a variety of different models [TT1 [19]. If we 
assume that this topologically disordered model is con- 
trolled by the same fixed point as the standard Ander- 
son model, then we expect Eq. [2] to hold in our system, 
with the same f(x). Previous work has expanded f(x) 
in polynomials in order to fit the data [HI [18] . We have 
not found this to be a successful strategy, at least in 
part because the underlying functions are strongly non- 
polynomial away from the transition region (see below). 
In general in such fits, there is a large number of pa- 
rameters [13 HH] , so it is helpful to constrain as many 
of them as possible. We assume universality and use 
four pieces of information from previous work: v = 1.58 
and /(0) = 0.576, as shown for the standard Ander- 
son model [17] . Further, we have asymptotic limits 
f(x —> ±00) = \x\ ±u [T5]. We choose to fit to logA lu = 
log/(0u> 1 ^) = h (4>w l l"). Then we have h (0) = -0.55 
and ho(x — >■ ±00) = ±vlog(±x). The asymptotic lim- 
its are obeyed by the function a(x) — i/sinh _1 (a;/2). 
We then fit our data to h (x) — a(x) + S(x) where 
(5(0) = -0.55 and S(\x\ -> 00) = 0. 

The scaling form, Eq. [2] applies only for w sufficiently 
large. At small w, corrections to scaling modify Eq. [2] 
and they have proven essential for accurate determina- 
tion of v and A c |17j . Corrections to scaling require intro- 
duction of some number of irrelevant variables r]i{E,p) 1 
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which have no effect on the scaling in the limit w — > oo. 
We illustrate with only one irrelevant variable ry, but 
more are easily added. We rewrite our scaling equation 
as \ogk w = h{(f)w 1 / 1 ' ,r]W' y '), where y is another critical 
exponent. We expand h in powers of rjw"^, 

\ogk w = rTw- m ^h m {<l>w x l v ), (3) 

m— 

where ho is the limiting one-parameter scaling function. 
It is often sufficient to keep only m = 0, 1 in Eq. [3j [17] 
which we will do here. We then define log A correctod = 
log A™ --qw-^h^w 1 /"). 

We are interested in finding the localization length at 
many values of E, unlike the usual choice of E = 0. 
The Green functions in the determination of A^, must be 
evaluated separately for each choice of E, increasing the 
computational requirements. We choose 14 values of den- 
sity and 40 values of energy satisfying —2 < E/Vo < 0.5, 
which are focused around the relevant energies for the 
metal insulator transition near p 1 / 3, rj 0.26. We do not 
study E = because isolated sites always produce eigen- 
values with E — 0, causing divergences in the Green func- 
tions. We take Lq = 15a* B and L c = mxa.[Lo,w/2]. At 
each p, we take w from 22a^ in increments of 3a* B up to 
the largest system size we find computationally tractable. 
For the 14 densities, the maximum w/a* B studied are 
(121, 103, 88, 79, 70, 61, 55, 52, 46, 43, 40, 37, 34, 31), from 
lowest to highest density. Previous studies of mobil- 
ity edges and critical exponents have performed separate 
scaling at each E [53]. In this work, all of the data col- 
lapse onto a single scaling curve, see Fig. [I] showing uni- 
versality independent of energy. These calculations used 
approximately 24 CPU-years of computing time. 

We approximate S(x) as a sum of three Gaussians, 
fixed to have 6(0) — —0.55, with eight free parame- 
ters. We use a direct search algorithm as implemented 
in MATLAB's fminsearch function to find the least- 
squares optimal choice of these eight parameters. Within 
the search loop, for each proposed S(x), we find (f>(p,E) 
(and possibly r](p, E)) independently at each value of 
(p,E), giving several hundred additional parameters. It 
would be preferable to have a functional form for 4>(p, E), 
but we did not find a parametrization that permitted 
accurate fits. To ensure that these fits are not overde- 
termined, we include in our dataset only (p, E) with at 
least 4 values of w; the average number of values of w at 
each (p, E) is 13. In fits including corrections to scaling, 
we approximate h\ (x) as a second order polynomial with 
hi(0) — 1, which fixes the scale of rj. 

Fig. [lji shows the scaled data, with one correction to 
scaling. The upper branch shows the extended states 
(A — > oo as w — > oo) and the lower branch shows lo- 
calized states. The scaling is excellent, which justifies 
the use of v and A c from the standard Anderson model 
studies. The consistency of A c with previous work indi- 
cates that the lattice spacing in the standard Anderson 



model is equivalent to a* B in this model. The extended 
state curve does not approach its <f> — ¥ oo asymptote 
due to a lack of data in the computationally demand- 
ing high density regime. Fig. [Tja shows the same fit with 
a and the irrelevant corrections subtracted. The dashed 
lines show the constituent Gaussians in 6, which are not 
well-constrained in the fits. See supplementary informa- 
tion for other fits. Fig. [2] shows the resultant correlation 
length £(p,E) = 101". At each p, the correlation length 
falls smoothly as E moves away from the derealization 
transition, just as we expect. These results give confi- 
dence that the 4>(fi, E) are not simply arbitrary parame- 
ters that happen to produce good scaling fits but rather 
are determined by the underlying physics. We find, how- 
ever, that 4>(p, E) (and thus £(p, E)) is quantitatively de- 
termined only for localized states. 

To assess the quality of the scaling curves produced, we 
consider the \ 2 statistic. Previous high accuracy numer- 
ical works on the standard Anderson model have evalu- 
ated their fits with the quality of fit Q (also known as the 
p- value) from the \ 2 statistic and the number of fitting 
parameters Np [T7]. All fits shown here have Q = 0, 
indicating that statistical fluctuations of the numerical 
procedure alone are insufficient to explain the deviations 
of the data from the model. This is not, however, sur- 
prising. Our study has several thousand degrees of free- 
dom, d = N B — Np — 1, where N B is the number of 
data points, which gives it statistical power to detect rel- 
atively small deviations between the model and the data. 
In the true scaling function, S(x) is not actually a sum of 
three Gaussians, and we are sensitive to the deviations 
between the true S(x) and its model. We should be able 
to add more Gaussians to 5 to better approach the uni- 
versal function, but fitting becomes difficult. If we are 
not entirely interested in the exact shape of 6 (x), then 
this deviation is not a concern. If we limit our dataset to 
contain only (p, E) close to the derealization transition, 
as in previous studies of the standard Anderson model, 
we can obtain good values of Q; this occurs because the 
number of data points is significantly reduced and the 
scaling function h(x) can be accurately approximated as 
a polynomial. 

Since we are interested in extracting correlation 
lengths £ = \<fi\ v , a better determination of the qual- 
ity and confidence of the fits is to compare the resul- 
tant 4>(p, E) for different fitting procedures. We compare 
fits with and without corrections to scaling and with the 
smallest w being 22, 25, or 28, in units of a* B . Depending 
on initial guesses, fits can arrive at a number of different 
local minima. We use a multistart procedure for the fit- 
ting, starting with 100 widely varying parameters for 5(x) 
and h\{x). We find that the "best fits," judged solely by 
minimizing x 2 , have highly oscillatory S(x), discontinu- 
ous (/>(p,E), large deviations from the asymptotic form 
even for strongly localized states, or large corrections to 
scaling; they generally have multiple of these features. 
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FIG. 2. (a) Correlation lengths £(p, £) from the fit of Fig. [T] 
for each of the 13 densities which produced enough data to be 
studied, offset for clarity. Solid circles mark localized states 
and crosses mark extended states. Solid and dashed lines 
are guides to the eye. As expected, the correlation length 
increases smoothly as the mobility edge is approached from 
either side. This figure shows the choice of energies for study, 
which are focused on the area of interest for the critical den- 
sity, near p 1 " = 0.24. The mobility edge is asymmetric just 
as is the density of states, with the upper mobility edge closer 
to E — than the lower mobility edge. The values of (p,E) in 
the lower-left corner of the figure have extremely low density 
of states, contributing less than 0.1% of total states, so are ex- 
cluded, (b) Deviations in fitted values of the scaling variable 
4>{p, E) from 92 different fits with and without corrections to 
scaling and with the smallest value of w taken to be 22, 25, 
or 28 (in units of a* B ). Localized states with cj> < —0.05 are 
quantitatively determined within 10% by the scaling fits. See 
the Supplementary Information for more versions of this plot. 



If we exclude the fits with these four characteristics, for 
w > 25 the best fits without (with) corrections to scaling 
have x 2 /g? ~ 70 (8). Including the anomalous fits, we 
can find x 2 /d k. 31 (4). We find cf>(p,E) from 92 differ- 
ent fits and, independently at each (p,E), find the mean 
(<ft) and standard deviation A(j>, shown in Fig. ^p. The 
values of <p(p, E) are found to vary by less than 10% in 
the localized regime, except for the points closest to the 
derealization transition. For £ < lOOa^, we can consider 
the localization lengths to be given quantitatively by the 
scaling method. Due to the lack of strongly-extended 
states, 5{x > 0) is not well-determined, and the fits show 
a range of different shapes. It is then no surprise that 
4>(p, E) in the delocalized region varies widely. Accumu- 
lating more data in the extended regime should fix this 



problem. 

This scaling technique quantitatively gives the local- 
ization lengths at nearly any localized (p, E) we care to 
study. Application of this method to systems with on- 
site disorder, in addition, should shift the mobility edges 
"inwards" so that more states are localized. These local- 
ization lengths allow insight into the properties of a range 
of material systems, and in future work we will consider 
their effects on intermediate band photovoltaics. 
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APPENDIX: Scaling curves without correc- 
tions to scaling 

As discussed in the main text, several fitting procedures 
can be used to produce plots similar to those in Figs. 1 
and 2. We illustrate here that the choice of fitting proce- 
dure does not significantly change the localization lengths 
extracted. Fig. [3] shows fits equivalent to those in Fig. 
1, without including corrections to scaling. The scaling 
still looks very good, though Fig. [3}d shows that the data 
points do not fall on the scaling curve as precisely as in 
Fig. lb. The overall shape of the scaling curve deviates 
from that in Fig. 1 mostly in the extended region, <j> > 0, 
which is not well determined from the data. Addition- 
ally, though the overall curve is very similar in shape to 
that in Fig. lb, the underlying Gaussians are different, 
which illustrates the difficulty in converging the fits of 
8{x). 

The correlation lengths £ extracted from these fits are 
shown in Fig. [4^. They are clearly very similar to those in 
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FIG. 3. Scaling function for data with w > 25, including no 
corrections to scaling. The scaling is very similar to that in 
Fig. 1 and produces similar £ (see Fig. [4^i) . In (b), blue points 
are the same as in Fig. lb while green points are those of the 
new fit. Dashed lines show the three Gaussians from the new 
fit of 8, which are clearly different from those in Fig. 1. There 
are 5829 data points and 449 fit parameters. The normalized 
X 2 statistic per degree of freedom d is 71.2. 39% of the data 
points are within error of the fit. 



FIG. 4. (a) Correlation lengths $(p,E) from the fit of Fig. 
|3J for each of the 13 densities which produced enough data 
to be studied, offset for clarity. Plot is similar to Fig. 2a. 
Solid circles mark localized states and crosses mark extended 
states. Solid and dashed lines are guides to the eye. (b) Sim- 
ilar to Fig. 2b, deviations in the fitted values of 4>{p, E). Blue 
points show deviations for 65 fits with no corrections to scal- 
ing. Green points show deviations for 27 fits with corrections 
to scaling. Both include data sets with w > 22, 25, and 28. 



Fig. 2a, showing that the scaling fits robustly determine 
£(p,E), regardless of whether corrections to scaling are 
used. Fig. 2b shows the deviations in the fitted values 
of for 92 different fits with and without corrections to 
scaling. 



Fits including corrections to scaling have larger dis- 
agreements in (f>(p,E). Fig. [4J3 shows the deviations for 
65 fits without corrections to scaling (blue) and 27 fits 
with corrections to scaling (green). The fits without cor- 
rections to scaling have a smaller number of fit param- 
eters, which seems to constrain the fits more. In any of 
the cases studied, the states with £ < 100a£ are quantita- 
tively determined within 10% by the fits. The consistency 
of the fits without corrections to scaling argues for use of 
that method in future studies, despite the large values of 
X 2 produced by ignoring the corrections to scaling. 



